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FIRST SEMI-ANNUAL FIRST YEAR PROGRESS REPORT 

This report covers technical progress during the first six months of the first year of 
NASA SR&T contract “Modeling the Magnetic and Thermal Structure of Active Regions,” 
NASW-03008, between NASA and Science Applications International Corporation, and covers 
the period January 14, 2003 to July 13, 2003. Under this contract SAIC has conducted 
research into theoretical modeling of the properties of active regions using the MHD model. 

Publication in The Astrophysical Journal 

The publication “Coronal Mass Ejection: Initiation, Magnetic Helicity, and Flux Ropes. 
I. Boundary Motion-Driven Evolution,” by T. Amari, J. F. Luciani, J. J. Aly, Z. Mikic, and J. A. 
Linker, has appeared in The Astrophysical Journal, 585, 1073 (2003). 

Abstract. In this paper we study a class of three-dimensional magnetohydrodynamic 
model problems that may be useful to understand the role of twisted flux ropes in coronal mass 
ejections. We construct in a half-space a series of force-free bipolar configurations with 
different helicity contents and bring them into an evolution by imposing to their footpoints on 
the boundary slow motions converging toward the inversion line. For all the cases that have 
been computed, this process leads, after a phase of quasi-static evolution, to the formation of a 
twisted flux rope by a reconnection process and to the global disruption of the configuration. 
In contrast with the results of some previous studies, however, the rope is never in equilibrium. 
It thus appears that the presence of a rope in the pre-eruptive phase is not a necessary condition 
for the disruption but may be the product of the disruption itself. Moreover, the helicity keeps 
an almost constant value during the evolution, and the problem of the origin of the helicity 
content of an eruptive configuration appears to be that of the initial force-free state. In addition 
to these numerical simulations, we report some new relations for the time variations of the 
energy and the magnetic helicity and develop a simple analytical model in which the magnetic 
field evolution exhibits essential features quite similar to those observed during the quasi-static 
phase in the numerics. 

A reprint of this paper is included in Appendix A. 

Presentation at the Spring AGU/EGS Meeting, Nice, France 

Zoran Mikic presented a paper at the Spring AGU/EGS Meeting, which was held in Nice, 
France, April 7-1 1, 2003: 

• Simulation ofCMEs Originating in Active Regions 
Z. Mikic, J. A. Linker, P. Riley, and R. Lionello 

Abstract. Previously we have addressed the initiation of CMEs by large-scale 
changes in the solar photospheric magnetic flux. These CMEs had a global nature, 
since their size was comparable to the solar radius. These early investigations were 
primarily focused on the theoretical aspects of CME initiation, and the models were 
thus rather idealized. We will describe recent advances in our computational models 
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that have extended our capability to study localized CMEs, such as those that might 
originate from an active region. This capability will allow us to analyze actual CME 
events, and to compare our results in detail with CME observations. 

This paper is included in Appendix B. 

Presentation at the Solar Physics Division Meeting, Columbia, Maryland 

Y ung Mok presented a paper at the Solar Physics Division Meeting, which was held in 

Columbia, Maryland, June 16-20, 2003: 

• Parametric Dependence of Coronal Heating Mechanisms and Active-Region 

Emissions 

Y. Mok, R. Lionello, Z. Mikic, and J. A. Linker 

Abstract. The thermal structure of an active region depends on the mechanism 
that heats the coronal plasma. A number of coronal heating mechanisms have been 
proposed over the years. They have different parametric dependences on the magnetic 
field, plasma density, and possibly other variables. Different mechanisms result in 
different thermal structures, and therefore, different EUV and soft X-ray emissions 
from an active region. Hence, the comparison between the computed emissions based 
on these models and the observed emissions will help to discover the parametric 
dependences of the actual heating mechanism and put some restrictions on the 
theoretical models. We have developed a 3D thermo-magnetohydrodynamic code to 
compute the thermal structure of an active region. The emissions resulted from 
various heating models will be compared with the images obtained from SOHO and 
Yohkoh. 

This paper is included in Appendix C. 
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APPENDIX A 

“Coronal Mass Ejection: Initiation, Magnetic Helicity, and Flux Ropes. 
1. Boundary Motion-Driven Evolution” 

T. Amari, J. F. Luciani, J. J. Aly, Z. Mikic, and J. A. Linker 
The Astrophysical Journal , 585, 1073 (2003) 
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ABSTRACT 

In this paper we study a class of three-dimensional magnetohydrodynamic model problems that may be 
useful to understand the role of twisted flux ropes in coronal mass ejections. We construct in a half-space a 
series of force-free bipolar configurations with different helicity contents and bring them into an evolution by 
imposing to their footpoints on the boundary slow motions converging toward the inversion line. For all the 
cases that have been computed, this process leads, after a phase of quasi-static evolution, to the formation of 
a twisted flux rope by a reconnection process and to the global disruption of the configuration. In contrast 
with the results of some previous studies, however, the rope is never in equilibrium. It thus appears that the 
presence of a rope in the preemptive phase is not a necessary condition for the disruption but may be the 
product of the disruption itself. Moreover, the helicity keeps an almost constant value during the evolution, 
and the problem of the origin of the helicity content of an eruptive configuration appears to be that of the 
initial force-free state. In addition to these numerical simulations, we report some new relations for the time 
variations of the energy and the magnetic helicity and develop a simple analytical model in which the 
magnetic field evolution exhibits essential features quite similar to those observed during the quasi-static 
phase in the numerics. 

Subject headings: MHD — stars: coronae — stars: flare — stars: magnetic fields — 

Sun: coronal mass ejections (CMEs) — Sun: flares 
On-line material : color figures 


1. INTRODUCTION 

The problem of the nature of the triggering of coronal 
mass ejections (CMEs) has attracted considerable interest in 
the last few years. Priest & Forbes (2001 ) present an interest- 
ing review of the many observations constraining a pre- 
emptive configuration. They stressed, in particular, two 
aspects that will be at the center of this paper: (1) shear is 
seen to be present in the chromosphere by the observations 
in Ha of the fibrils (Tanaka & Nakagawa 1973) and along 
the inversion line of the normal component of the magnetic 
field (Hagyard et al. 1982; Machado & Moore 1991; 
Hagyard 1990), and (2) converging motions toward the 
inversion line have been observed, and they may increase 
the shear in a presheared magnetic configuration. 

Another important feature observed in a relatively large 
number of CMEs is the presence of a prominence (Priest & 
Forbes 2001 ) and the ejection of a plasmoid, probably at the 
origin of the magnetic cloud sometimes observed afterward 
in the interplanetary medium (Burgala et al. 1981; van 
Driel-Gesztelyi, Schmieder, & Baranyi 2002). These coronal 
structures contribute in different ways to the total magnetic 
helicity. Twisted flux ropes have been often thought to be 
good candidates for explaining the coronal helicity content 
and budget and for triggering the eruptive events (Low 
1994; Amari et al. 1999c, 2000). However, one important 

1 CEA/DSM/DAPNIA, Service d’Astrophysique (URA 2052 associee 
au CNRS), Centre d’Etudes de Saclay, F-91 191 Gif-sur- Yvette Cedex, 
France. 


issue has not been settled yet: is it really necessary to have a 
twisted flux rope (in equilibrium) prior to the disruption 
(Amari et al. 2000; Linker et al. 2001), or is the latter created 
as a consequence of reconnection during the disruption 
itself (Antiochos, Devore, & Klimchuk 1999). It is worth 
noticing that these considerations do not exclude a priori 
one or the other topology since both the sheared arcade and 
the twisted flux rope may be sources of magnetic helicity 
exchange between the convection zone and the corona 
(there is conservation of the total helicity contained in the 
domain convection zone + corona, which is filled up with 
highly conducting plasma). 

Modelizations of the triggering of CMEs have mostly 
concentrated on the problem of the evolution of a magnetic 
field occupying some unbounded domain (generally either a 
half-space or the exterior of a sphere) and embedded in very 
tenuous highly conducting plasma, this evolution being 
driven by slow (compared to the Alfven velocity) motions of 
the footpoints. Initially, the efforts have concentrated on 
axisymmetric arcades submitted to purely shearing boun- 
dary motions. This problem has been treated in both Carte- 
sian (Aly & Amari 1985; Aly 1990; Amari et al. 1996a, 1997 
and references therein) and spherical geometry (Mikic & 
Linker 1994; Aly 1995), the evolution being found in each 
case to lead to the formation and ejection of a plasmoid 
once a sufficient shear has been applied. The effects of con- 
verging motions have also been considered analytically in 
the Cartesian case, and the existence of a catastrophic non- 
equilibrium transition was obtained (Priest & Forbes 1990; 


1073 


1074 


AMARI ET AL. 


Vol. 585 


Forbes & Priest 1995). By running resistive simulations, it 
was also proved in Forbes (1991) and Inhester, Birn, & 
Hesse (1992) that a plasmoid gets formed and that the proc- 
ess exhibits an impulsive phase. There is, however, one 
strong limitation inherent to all these two-dimensional stud- 
ies: the produced flux rope is not anchored in the boundary. 
Then the necessity of performing three-dimensional calcula- 
tions was strongly underlined (Priest & Forbes 2001 ). 

Three-dimensional approaches have been mainly con- 
cerned with the effects of a shear applied to the footpoints of 
an initial bipolar magnetic configuration. To avoid flux 
pileup on the walls of the computational box, it proved nec- 
essary to give a twisting component to these motions ( Amari 
et al. 1996b; Tokman & Bellan 2002). However, this had the 
bad effect of limiting the amount of shear along the neutral 
line, which turned out to be not as strong and coherent as in 
the axisymmetric case (Hagyard 1990). Consequently, it 
was noted (Amari et al. 1996b) that some of the field lines 
may lean sideway, letting the inner flux rope merge through 
them and limiting below the amount of current associated 
with the current sheet that would form in two dimensions. A 
transition toward partial opening was obtained, but without 
strong dissipation and ejection of a plasmoid. 

The effect of converging boundary motions has not been 
considered yet, and it is the aim of this paper to investigate 
their possible effects on a bipolar configuration (complex 
topology fields will not be introduced here) initially in a 
force-free state and having nonzero helicity contents. In the 
solar corona, such configurations may be the result of sev- 
eral different processes. They may be associated with a 
newly emerging active region bringing electric current with 
it (Leka et al. 1996), or they may be the remnant of poster- 
uptive arcades that have not relaxed to a potential state. 
(According to Taylor’s classical conjecture, it could be sup- 
posed that the field relaxes to a constant-a force-free config- 
uration after an eruptive event. Recently, however, the 
validity of this idea has been seriously challenged. It has 
been found indeed to be incorrect in numerical MHD simu- 
lations of coronal disruption [Amari & Luciani 2000] and in 
an analysis of the global budget of an observed flaring active 
region ejecting at infinity in a magnetic cloud some amount 
of helicity [Bleybel et al. 2002]). However, the details of 
these processes, as well as the mechanisms that could lead to 
the formation in the magnetic structure of a prominence by 
condensation of cold material, are of no concern to us here. 

For our purpose, we first construct numerically in the half- 
space Cl = {z > 0} (or more exactly in a cubic box of very 
large size) a series of initial force-free configurations with the 
adequate structure. Next, we impose on the boundary 
{z = 0} some motions converging toward the inversion line, 
and we determine the resulting evolution of each of our initial 
configurations by solving the full set of MHD equations. 
Some of the important questions that we may answer that 
way are as follows: ( 1 ) Do the converging motions contribute 
to the helicity contents of the magnetic structure? (2) How 
long can the field evolve quietly in an approximately quasi- 
static way? (3) What happens when the quiet phase ends, is 
there production of a twisted flux rope in equilibrium, or is 
the system subject immediately to a global disruption? 

The paper is organized as follows: First, we present some 
general analytical results. We recall in § 2 (after a descrip- 
tion of our model) some properties of two basic physical 
quantities, the magnetic energy and the relative magnetic 
helicity, and we derive in § 3 some new formulae governing 


their evolution in a quasi-ideal situation. We consider in § 4 
a simple exact solution describing, in the quasi-static frame- 
work, the evolution of a force-free field driven by a special 
class of boundary converging motions. The second and 
main part of the paper is devoted to a description of the 
results of our numerical simulations. The initial force-free 
states are constructed in § 5, and their evolution driven by 
converging boundary motions is considered in § 6. Finally, 
we discuss the significance of our results, in particular their 
relevance to a theory of CME, in the concluding § 7. 

2. MODEL PROBLEM: A CLASS OF INITIAL 
BOUNDARY VALUE PROBLEMS 

This section is devoted to a short description of our model 
and to a reminder of some basic properties of the magnetic 
energy and the relative magnetic helicity. 

2. 1 . Description of the Model 

In the model we consider in this paper, the corona and the 
photosphere are represented by the half-space O = {z > 0} 
and the plane S = {z = 0}, respectively. Cl is assumed to be 
filled up with a low-/? slightly resistive and viscous plasma 
embedded in a magnetic field B(r , /). At some initial time to , 
the field B () {r) = B(r, to) is taken to be force-free and to have 
a finite energy W(to) and a finite magnetic relative helicity 
H(to). Therefore, it does obey the equations 

VXBo = a 0 Bo , (1) 

V • Z?o = 0 , (2) 

with ao(r) satisfying the well-known constraint 

Bo • Va 0 = 0 . (3) 

In practice, this initial force-free field can be obtained by 
several different procedures. Hereafter we shall use a numer- 
ical twist and relaxation method, but we could also con- 
struct Bo by solving a boundary value problem (BVP) in 
which the normal component B z and the function a are 
specified, respectively, on the whole boundary 5, and on 
that part S + of S where B z > 0 as in Sakurai (1981) and 
Amari, Boulmezaoud, & Mikic (1999a). Although no math- 
ematical existence theorem has yet been proven for this 
problem, it seems quite likely that it has always a solution, 
at least if a is chosen not too large. Indeed, the existence of 
solutions has been now established for the similar problem 
set: either ( 1 ) in a bounded domain (Boulmezaoud & Amari 
2000), this result being directly relevant in the situation 
actually considered in our numerics, which are done not in 
Cl, but in a cube Cl/,; or (2) in an unbounded one of the “ exte- 
rior type,” like the exterior of a sphere (Kaiser, Neudert, & 
von Wahl 2000). 

This initial configuration being given, we start imposing 
for / > /o some slow motions to its footpoints on the boun- 
dary S, and the field is thus driven into an evolution. These 
motions are parallel to the surface, i.e., plasma is not 
allowed to go through 5, and the given velocity v(x,y, 0, t) 
thus satisfies 

V;(x,y,0,t) = 0 . (4) 

Of course, our aim is to determine the evolution of the sys- 
tern, i.e., to compute the field B(r , /), and we shall do that by 
solving the full set of equations of magnetohydrodynamics. 
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2.2. Magnetic Energy 


and 


The magnetic energy of the field at time t is defined by 

w{t)= hj i / [t)dv - (5) 

This quantity bears some interesting relationships with 
the energies of two auxiliary fields that play an important 
role in all the evolutionary problems of the type considered 
in this paper: the potential field B n (t) and the open field 
B a (t), which are uniquely defined by imposing at time t the 
boundary conditions 


Therefore, A n: is a harmonic function that satisfies a homo- 
geneous Neumann boundary condition on S and vanishes 
at infinity, which implies 

A nz = 0 in ft . (18) 

Moreover, we can deduce from equations (14) and (15) that 
there does exist a potential x(x, y, 0 such that 


BK Z (x,y,0,t) = B az (x,y,0,t) = B : (x,y,0,t) (6) 


A s = A ns = V^X z on S (19) 


and an asymptotic decay at infinity. These fields have ener- 
gies given, respectively, by the standard relations 

16 t t-Jsxs w-n 

w a (t) = JL f \ B y^yMB : {Ay’At)\ dsdsl (8) 

l6ir 2 Jsxs k-r' | 

As it is well known, the energy W(t) always satisfies 

fVAO < mo • (9) 


Moreover, it is expected on the basis of quite general argu- 
ments (Aly 1984, 1991; Sturrock 1991) that 


W(t 0 )< W a (t 0 ) , (10) 

this inequality staying actually valid for t > t 0 as long as the 
field B(t) approximately evolves in a quasi-static way 
through a sequence of force-free configurations. 


and 


B z = -Vf x on 5 . (20) 

With the gauges fixed as above, the relative helicity is given 
by 


H = f (A • B - A* • B n )dv . (21) 

Jn 

Before closing this subsection, we add one remark of 
mathematical character. We would like to stress that a 
rigorous theory of the vector potential in the half-space ft 
has been recently worked out in Boulmezaoud (1999). 
Boulmezaoud defines for any real 7 the weighted space [with 
weight function p(r)\ 


H 1 (div; ft) = 
&<SI) 


v e 


-\ 3 


P~ l (r)v_ 


€ [L 2 (n)] 3 ,^ +1 (r)V.»G [L 2 (n)] 3 [ , (22) 


2.3. Magnetic Helicity 

In the calculations below, we shall follow in particular the 
evolution of the relative magnetic helicity H(t). Then we 
recall here the definition of that important quantity (Berger 
& Field 1984). We first introduce the representations 

B — \ X A (11) 

and 

B n = \XA n = \Vn (12) 

of B and B n in terms of the potential vectors A(t) and An{t) 
and the scalar potential V n (t). 

Using the gauge arbitrariness of A and A n , we can choose 
these quantities in such a way that, at any time t y 

V-A = V-An = 0 in ft, (13) 

A s = A^ on S, (14) 

\ s 'A s = V* •Ans = 0 on 5, (15) 

where X s denotes quite generally the component of 
X(x,y, 0, t) parallel to the boundary S and = d x x -f d y y. 
The compatibility of these conditions is easily verified. We 
just note that imposing equation (14) is possible because of 
equation (6). From equations ( 1 3) and ( 1 5), we get 

-■ V 2 A„ = [v X (V X A r ) - V(V . A „)) . z 

= (VXfi I ).z = 0 (16) 


which is a Hilbert space for the norm 

IMk(div;«)= [l|P 7 ('-)r|lo,n+||P 7+, WV.« > ||y V2 , (23) 

and shows the existence of A e [ (ft )] 3 (in an integer) 
satisfying \ • A = 0 in ft, if either B e (ft)] 3 , with 
m > 0, or B € H\ (div; ft), with m > 1 . 

3. EVOLUTION OF THE MAGNETIC HELICITY AND 
THE MAGNETIC ENERGY 

We now establish some general formulae for the time var- 
iations of the magnetic energy and the magnetic helicity, 
assuming that we are in a phase of evolution where dissipa- 
tive processes (either resistive or viscous) can be neglected. 

3.1. The Tangential Electric Field on S 

With resistive effects negligible on S and boundary 
motions satisfying equation (4), we have from Ohm’s law 

—cz X E s = zX ( vX B) s = B z v s on 5 . (24) 

By a well-known theorem, the quantity on the right-hand 
side of the latter relation can be decomposed into an irrota- 
tional part and a solenoidal one (Helmholtz’s decomposi- 
tion of a two-dimensional vector). More precisely, there do 
exist two functions/(x,^, /) and g(x,y , t) such that 

B z v s = V/ + \ s gX z • 


(25) 
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Here /and g need to satisfy the equations 

V 2 J = V s .(B : v s ), (26) 

V^ = -[V,.X (B : v t )]-z, (27) 

obtained by taking successively the surface divergence and 
rotational of equation (25), and they are uniquely deter- 
mined if we impose in addition that they do decay to 0 at 
infinity on S . 

Using Faraday’s induction law, we also have 

^ = -cV • (E s X z) = -V . (B : v s ) = -V]f , (28) 

and we obtain by comparing this relation with equation (20) 

d\ 

< 29 > 

3.2. Magnetic Energy 

The evolution of the magnetic energy of the configuration 
is given by the flux of the Poynting vector through the boun- 
dary S, i.e., 

W = ^-J{E s XB s ).zds. (30) 

This relation can be given several equivalent forms by using 
the relations of the previous subsection. First, we get with 
the help of equation (24) 

W = -±-j(v s -B s )B z d s. (31) 

Next, we use the decomposition given by equation (25) in 
the latter relation, which leads to 

w = - T j B . (V/ + \,g X Z)ds 

= “///• <JB, + gzx b s ) -f\ s • b s 
+ gz • x B s ]ds 

< J2 > 

Here we have effected an integration by parts and used 
Gauss’s theorem on S, equation (2), and Ampere’s law 

VXB = ^j, (33) 

where j is the electric current density. In equation (32), W 
appears to be divided into two terms: one related to / and 
then to the change of B z . , and another one related to g and 
then to the rotational part of B z v s . 

In the case in which the boundary motions keep B z invari- 
ant if = 0) and the field in Q is force-free (4 nj z = acB z ), 
equation (32) reduces to 

w = -^- [gaB.ds, (34) 

47 rJs 

a formula that was derived for the first time in Aly (1991). 


3.3. Magnetic Helicity 

The variation of the relative helicity in an arbitrary 
domain D is given by (Berger & Field 1984) 

H = 2 [ ( cA„ X E s - d, V^A*) *hds , (35) 

JdD 

where h is the exterior unit normal to the boundary dD of D 
and the potential vectors of B and B- are assumed to satisfy 
the gauge constraints given by equations (13) and (14). With 
D = Q and the gauge condition given by equation (15) 
enforced, we have A rz = 0 by equation (18), and the for- 
mula above reduces to 

H = -2 j A ws • (cE s X z)ds . (36) 

Using the representation given by equation (19) for A ns 
and equation (24) for E s in the latter equation leads to 

H = -2 J(\ s xX z) • (B~v s )ds , (37) 

whence, by injecting the decomposition given by equation 
(25), 

H = - 2 f s ^ sX X + x *) ds 

= -2 / (VJX\ sX )-zds-2 [ y s g • ds 

Js Js 

= -2 JjV s *(fVsX)}-zds 
-2 h.-0fl.x)* + ij^Vlxds. ( 38 ) 

The first and second integrals on the right-hand side disap- 
pear by Stokes and Gauss’s theorems, respectively, and the 
assumed asymptotic decay of the various functions. There- 
fore, by taking equation (20) into account, we are finally left 
with the simple expression 

H = -2 J gB : ds . (39) 

Of the two functions /and g , only g (which is determined by 
the rotational part of B z v s , which does not change B z ) 
appears explicitly in this formula. This does not mean, how- 
ever, that / does not play any role: actually, its effect is con- 
tained in the time dependence of B z given by equation (28). 

Another form of equation (39) is obtained by intro- 
ducing the flow generated by the velocity field v s on S. If we 
denote as r = r(i*o, t) the position at time t of a material 
point located at ro at the initial time to, we have 
B z (r,t)ds(r,t) = B 0z (ro)dso(ro) by flux conservation and 
then 

H = ~2 J g[r(r 0 ,t),t]Bo : (r 0 )dso . (40) 

In this form, all the time dependence has been transferred to 
the function g, now expressed in terms of Lagrangian coor- 
dinates rather than Eulerian ones. 

It is worth emphasizing here an important point. The 
change in the magnetic helicity depends only on quantities 
(B : and v on S) that are given in our problem as boundary 


No. 2, 2003 


CORONAL MASS EJECTION INITIATION. I. 


1077 


conditions and not at all on the dynamics of the field B 
inside ft. H(t) can thus be computed a priori, without having 
to solve the MHD equations. This makes a striking differ- 
ence with the energy, whose variation (e.g., eq. [30]) depends 
on the tangential component of B on 5, a quantity that has 
to be computed a posteriori from the solution B(r , t) in ft. 

3.4. Symmetric Situations 

Let us now consider how the helicity evolves when the 
boundary conditions satisfy some symmetry conditions. As 
a first example, let us assume that, like in our numerics 
reported hereafter, 

B z (x, —y, 0, to) = -B z (x,y,0,t 0 ) , (41) 

and either 

»*(*> 0, t) = - v x (x,y , 0, t) , (42) 

v y (x,-y,to,t) = v y (x,y,0,t) , (43) 

or 

v x (x,-y,0,t) = v x (x,y,0,t) , (44) 

v y (x,-y,0,t) = -v y (x,y,0,t) . (45) 


4. A SIMPLE ANALYTICAL MODEL 

In this section we present a simple analytical model 
describing the quasi-static evolution of a force-free mag- 
netic field whose footpoints are submitted to a particular 
class of converging or diverging motions. Comparison of 
the behaviors of the fields obtained in this model and in the 
dynamical calculations reported in the next section will lead 
to some interesting conclusions. 

4 A. A Sequence of Force-free Fields 

Let us choose some continuously differentiable positive 
function A(/) satisfying A(/o) = 1 and define for t > to the 
time sequence of fields 

B(r,t) = \ 2 (t)B 0 [\m, (51) 

where Bo is the initial force-free field introduced in § 2.1. 
The fields B(t) are also defined in ft, where they are easily 
checked to satisfy 

VB= 0, (52) 

\XB=aB, (53) 

with 


In the first case, we impose in addition that / = 0. Then 
B- is time independent on S , and equation (41 ) is also satis- 
fied at t > to . Using equation (27), we get 

(gB : ){x, -y) = ( gB : )(x,y ) (46) 

and 

H = -4 J gB : ds (47) 

(B-_> 0 on S + c S'). Therefore, we obtain a monotonically 
increasing (decreasing) helicity if we choose g to be negative 
(positive) on S + . Equation (47) does apply with a v of the 
vortex type considered in § 5.3. However, it is not valid 
when v is of the purely shearing type v = v(y)x, with 
v(-y) = -v{y). Such a v [used, e.g., in Devore 2000, with 
y(y) = uy] leads indeed to a changing B. [clearly, 
B s (x,y, 0, t) = B z {x,y - v(y)(t - t 0 ), 0, / 0 )]. 

In the second case (which is that of the converging 
motions defined in § 3.1), equation (41) is automatically 
conserved in time. We have 


a{r, t) = A(/)a o [A(0»-] . (54) 

Then all the fields of the sequence turn out to be force-free 
magnetic fields. Moreover, if we set 

< 55) 

B satisfies the equation of evolution 
r) R 

— = VX(vXB) (56) 

valid for a field embedded in a perfectly conducting plasma. 
Therefore, as the velocity v(x,y, 0, t) is parallel to 5, it is pos- 
sible to interpret the sequence as describing a quasi-static 
evolution of B driven by footpoint motions. In particular, 
we see that imposing A(0 to be monotonically increasing 
(decreasing) allows us to describe an evolution driven by 
motions converging toward (diverging from) the origin O. 

4.2. Some Properties of the Sequence 


(gB : ){x, —y) = —(gB : )(x,y) 


and 


(48) 


On 5, the normal component of B is given by 

**(*,*0,0 = A 2 (/)Z?o_-[A(/)x, A(0*0] , (57) 


H = 0 , (49) 

i.e., the helicity keeps a constant value. 

As a second example, let us consider the case of an axi- 
symmetric field. Then f g , and B z on S depend only on r, 
and rB z — A', with A the standard flux function. Whence 

H = 2 J ^g'ds = -2 J jAB : ds , (50) 

where we have effected an integration by parts with respect 
to r and used equation (25) to evaluate g'(r). Then we find 
(up to a different normalization factor) a relation derived in 
Aly (1992; his eq. [1 16]) under the supplementary assump- 
tion (not made here) v r = 0 on S. 


and because of equation (20), we have 

= x[A(<)Mo] = Xo[A(/)r] . 
Using equation (29), we thus get 

fir (*•> 0 = M‘)r ‘ ( V iXo)[A(/)r) 
A, A 

~ A (/) r dr r ’ * ’ 
or 


/= 


A 

A 


r 


dx 

dr ’ 


(58) 


(59) 


(60) 
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Moreover, we have from equation (25) 


v/ + V^x z = -jB : r , 

whose orthoradial component writes 

!^_% = 0 

r d(j) dr 

Along with equation (60), this gives 


( 61 ) 

(62) 


Xdx 

A <90 ’ 


(63) 


Both / and g can thus be expressed in terms of A and x, or xo 
by using equation (58). Note that we also have from 
equation (59) 


from ro to r(ro) [by admissible, we mean here one-to-one, of 
Jacobian J = D(r)/D(ro) > 0, and not moving the boun- 
dary points at which B z ± 0] and by comparing the energy 

W = ^j^\ B -\ a rUv (74) 

of the deformed field with that of the equilibrium field. 
Noticing that we may establish a one-to-one correspond- 
ence between the deformations i*(i*o) admissible for Bq and 
those r,(ro) admissible for B(t) by setting 

r(r 0 )< — >r,(r 0 ) = jLr[A(r> 0 ] , (75) 

we get for two deformations associated that way 


dx 

dt 



VsX 


0 , 


(64) 


i.e., x is a Lagrangian invariant on S. 

As an immediate consequence of equation (6), the poten- 
tial field B z (t) and the open field B a (t ), which are fully deter- 
mined by this boundary condition, are given by 

=A 2 (0/MA(f)r,0], (65) 

B„{r, t) =A 2 (/)£„[A(f)r,0]. (66) 


Making a simple change of variable, it is seen at once that 
the respective energies of B(t), B n (t), and B a (t) vary accord- 
ing to 


W(t) =X(t)lV(t 0 ), (67) 

fVAO = A(r) Wjr(to) , (68) 

W a (t) = \(t)W a {t 0 ) . (69) 

All these energies increase (decrease) monotonically if the 
motions are converging (diverging). Moreover, we have 

W(t) _ W^t) W a {t) 

W(*o) W n (t 0 ) W a {t 0 ) ’ 1 ' 

which implies that the energy of the open field cannot be 
exceeded at some time t if it is not exceeded at the initial time 
tQ. Then the sequence would be useless for someone wanting 
to disprove the validity of the conjecture quoted in § 2.2, 
according to which one has always W(t) < W a (t) when the 
field is force-free. 

Consider now the relative helicity of the fields B(t). Mak- 
ing the obviously possible choices 

A(r,t) =A(MA(r)r,/ 0 ], (71) 

4 T (r,t) = A(/)/4„[A(/)r, /o] , (72) 

for the potential vectors of B(t) and B~(t), respectively, we 
obtain 


W(t) - W{t) = A(/)[H^(0) - W(0)] , (76) 

which shows clearly that B(t) has the same stability proper- 
ties as Z?(0). If the latter is stable (unstable), so is /?(/). Note 
that we have not supposed the displacements |r(r 0 ) — fo| to 
be arbitrarily small, and then stability may be understood in 
the statement above as being either linear or nonlinear. In 
particular, we can claim the following result. Suppose that 
we start from a force-free field Z?(0) that is unconditionally 
nonlinearly stable (i.e., it is an absolute minimum energy 
state among all the configurations that can be obtained 
from it by an ideal finite deformation): then B(t) has the 
same property. 

As far as the specific problem considered in this paper is 
concerned, we can essentially retain the following statement 
from the considerations above: from an initial force-free 
state B 0 and a monotonically increasing function A(/) satis- 
fying A(/o) = 1, it is possible to construct a time sequence of 
force-free configurations B(t) such that (1) B(t) describes an 
ideal quasi-static evolution starting from B(t 0 ) = Bo and 
driven by converging boundary motions on S, (2) all the 
fields B(t) are as (ideally) stable as Bo, (3) the relative helic- 
ity of B(t) is equal to that of Bo, and (4) the energy of B(t) 
increases monotonically in time. 


5. CONSTRUCTION OF A NONZERO RELATIVE 
HELICITY FORCE-FREE CONFIGURATION B () 

We now start presenting our numerical simulations, 
which are done by substituting to the half-space 0 a cubic 
box of large size (compared to the characteristic scale of 
our boundary conditions). As a first step, we construct in 
this section a series of initial force-free states Bo. For that, 
we set up at time / = 0 a potential field, which we submit to 
twisting boundary motions for t e [0, t s ] and to a viscous 
relaxation, with no boundary motions, for / G [/*., / 0 ]. 


H{t) = 7/(0) (73) 

(just effect a change of variable in the integral defining H). 
Then the helicity keeps a constant value all along the 
sequence. 

Finally, let us consider the ideal stability properties of the 
fields B(t). We may test the stability of an equilibrium by 
applying to it an admissible displacement moving a point 


5. 1 . Initial Potential Configuration 
Our initial potential magnetic field 

B n (x,y, z, 0) = \ V n (x,y, z, 0) (77) 

is computed by solving the following mixed (Neumann- 
Dirichlet) BVP for the scalar potential V* G ^"^(Q/,), 
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V 2 K ff = 0 in n h , 

(78) 

= for z = 0 , 

(79) 

K = C on dtt/,\{z = 0} , 

(80) 


where £ represents the given distribution of B- on {: = 0} 
and 


CM = - 


J_ [ CM) 

2tr Js |r-r'| 


c/s' 


(81) 


is the value that would be taken by V n on the boundary of 
the box (z ^ 0) if the calculation was done in the whole Cl. 
Here we take 


N 

&x,y) = ]T a+e r ^ xy) + a k .e r ^ M , (82) 

A'=l 

with 

(83) 

a x °y 

Although this distribution of flux does not have a compact 
support, it is clearly decreasing fast enough in U 2 for equa- 
tion (81) to be applied. As for the values of the parameters, 
we choose = [—20, 20] x [—20, 20] x [0, 40) for our 

computational domain and set N — 1 , of = 2; erf = 1 , 

a t = ~ a k 7 *; and <i = x 7i = °. yti = -yj\ = °- 8 - As 
shown in Figure 1, the field Z?~ constructed that way mimics 
an arcade-like structure in the upper half-space ft. Note that 
B~ represents a good numerical equilibrium (electric current 
and magnetic forces are numerically small), since it has been 
constructed as the solution of a BVP on our computational 
mesh rather than by effecting a viscous relaxation as in 
Amari et al. (1996b). 



Fig. 1.— Selected field lines of the initial bipolar potential magnetic con- 
figuration. [See the electronic edition of the Journal for a color version of this 
fhjure.) 


5.2. M HD Equations 

Before describing the deformation of into a nonlinear 
force-free state at / o> let us write the system of the MHD 
equations we shall use. In nondimensionalized form, it is 
given by (Amari etal. 1996b; Amari, Luciani, & Joly 1999b) 

P^= -p(».Vc) + (VX B)XB-Xp 


+\ • (vpXv) + pg . (84) 

f)B 

— = XX(vXB)-XX(r,j), (85) 

§7 •(/»), ( 86 ) 

^=-„.Vp-rp(V. v) + H, (87) 

j = \XB (88) 

V . B = 0 , (89) 


where p, u, 77, and T are the mass density, the kinematic vis- 
cosity, the resistivity, and the adiabatic index of the plasma, 
respectively. 

Small values are used in this paper for the dissipation 
coefficients: v = 10 -2 to 10“ 3 for the kinematic viscosity and 
77 = 10~ 4 , 10~ 5 , 0 for the resistivity (for our mesh resolution, 
this gives Lundquist numbers of order 10 4 , 10 5 ), which 
allows us to neglect the term H in equation (87). The plasma 
P is taken to be 10~ 3 (i.e., of the order of the very small 
values observed in the corona) or 0, without any differences 
being actually found between the results. When we choose 
(3 — 0, we need to fix arbitrarily a mass density profile and 
of course neglect the gravity term in equation (84). Here we 
choose p = B 2 , which ensures a constant Alfven velocity, or 
p = 1 . Alternative choices of density profiles (exhibiting, for 
instance, a slower decrease with distance) have been found 
not to lead to sensibly different results. 

The MHD equations are discretized on a nonuniform 
mesh (141 x 121 x 91 nodes) and solved by using our semi- 
implicit scheme (Amari et al. 1999b). 

5.3. Twisting Evolution: f = 0 

For t > 0, we impose a twisting velocity field to the foot- 
points at the bottom of the box. The velocity field is pre- 
scribed by choosing in equation (25) 


/ = 0 (90) 

(therefore the distribution of B z is preserved for z = 0), and 
to allow comparison with our previous calculations, we use 
as in Amari et al. (1996b) the velocity field v s = V^X z, 
with 

(j> = <t>(B : , t) = D 0 R(/)Bre _ K fl ^ B ^“) /B:m “] , (91 ) 

where B :m ^ = sup s £ and vo = 10 -2 (then tfo is small com- 
pared to the Alfven speed v A = 1). R is a ramp function that 
is used to smoothly switch on or off the velocity field. It is 
chosen to be linear, the maximum velocity 10 -2 being 
reached at / = 10r A . This velocity field corresponds to two 
parallel vortices rotating in the same direction, as in Amari 
et al. (1996b). It is shown in Figure 2. Although an explicit 
formula for g is not given, it appears clearly that one may 
solve equation (27) for the unique solution p, which we will 
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— > ■ 1 .000e+00 (a) (b) 

Fig. 2. — {a) Vector plot and ( b ) cut at {x = 0} of the twisting boundary flow applied to the footpoints of the field lines of the initial configuration shown in 
Fig. 1 used to construct a nonzero helicity force-free field. 


nevertheless not need to use explicitly in what follows. Since 
from the definition of g one may deduce g'(B z ) = B z (t>'(B z ), 
one may easily check that the parity conditions of equation 
(46) for g are fulfilled. On the other faces of the box, we take 
the homogeneous Dirichlet condition v = 0 as in Amari et 
al. (1996b), which is the adequate condition for a viscous 
plasma in contact with a nonmoving wall. The viscosity is 
fixed to 10" 3 and the uniform resistivity vanishes (rj = 0). 

The photospheric twisting motions are applied from 
t = 0 to / = t s (with these times being expressed in units 
of t a ). By giving to t s different values (here 
t s = {0, 50, 100,200,400}), it is possible to generate a series 
of states U ts = (B , *,vt*,pt*,p t *) ranging from moderately 
sheared to highly sheared configurations. We use U° as a 
reference state. As in Amari et al. (1996b), the shearing 
process brings up shear near the neutral line while twist is 
introduced as one moves away from the neutral line on 
{z = 0}. The evolution along this phase is found to be 
almost quasi-static. The magnetic field acquires an energy 
that is monotonically increasing (see eq. [32] with / = 0) 
and magnetic helicity. Consistently with equation (47) (we 
are here in a situation where eqs. [41], [42], and [43] are satis- 
fied), the helicity turns out to be negative, its absolute value 
increasing monotonically and even at a constant rate after 
the transient phase during which the ramp is applied. As 
shown in Figures 3 and 4, \H^\ and W u are increasing 
functions of t s . 

5.4. Relaxation to an Equilibrium: f = g = 0 

At t = t s the driving photospheric velocity is switched off 
and the system evolves ideally (rj = 0) through a viscous 
relaxation procedure that is stopped at t = /o, when the field 
has reached an accessible stable numerical force-free state. 

As expected from equations (32) and (39), this phase 
(f = g = 0) implies no injection of magnetic energy and 


W0/Wrc 

W050/Wti 

WIOO/Wt: 

W200/W n 

W400/Wrr 



Fig. 3. — Variation of the magnetic energy for various periods during 
which the twisting boundary flow on Fig. 2 is applied to the potential con- 
figuration from / = 0 to / = t s and then switched off to reach a relaxed state 
at t = to, and then converging motions are applied. The various cases 
t s = 0, 50, 100, 200, 400 (in units of r A ) are represented. Magnetic energy W 
is expressed in units of W r (see text). The arrows indicate the three phases 
for the particular case t s = 400 in the same units. [See the electronic edition 
of the Journal for a color version of this figure .] 
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H100 

H200 

H400 



Fig. 4. — Variation of the magnetic helicity for various periods during 
which the twisting boundary flow on Fig. 2 is applied to the potential con- 
figuration from 0 to t — t s and then switched off for reaching a relaxed state 
at i — to, and then converging motions are applied. The various cases 
h — 0. 50, 100, 200, 400 (in units of r A ) are represented. Magnetic helicity 
H represents the relative helicity with our gauge (see text). The arrow indi- 
cates the starting beyond which converging motions are applied for the 
particular case t s = 200. [See the electronic edition of the Journal for a color 
version of this figure.] 


helicity. That this is the case indeed is clearly seen in Figures 
3 and 4, where these quantities appear to keep constant val- 
ues for t e [/ v , Jo] and for all the values of t s , apart from a 
small variation associated with the linear ramp. Various 
equilibrium states U ts with various energy and magnetic hel- 
icity levels are thus obtained at / = Jo- 
lt is worth noticing that we have chosen, on purpose, the 
times t s small enough in order to not reach the critical 
threshold beyond which the configuration starts to inflate 
more rapidly (Amari et al. 1996b). Otherwise, we would not 
get equilibria at / 0 . Figure 5 shows the equilibrium configu- 
rations obtained at l s = 200t a and to = 400t a . We have 
selected three subsets of field lines to illustrate the effect of 
the applied velocity field (shear and rotation). The amount 
of magnetic helicity injected in the initial configuration is 
clearly not small. We remark in this example that the mag- 
netic energy, although not negligible compared to the 
potential energy, is still below the energy of the open field 
which we found to be W a (t) / W n (0) = 3 for t < to 
during this phase that preserves B z , upon which W n 
depends. 


6. CONVERGING MOTION-DRIVEN 
EVOLUTION: / ^ 0 

We now apply for t > to converging motions to the foot- 
points of the force-free fields = B(to) constructed in the 
previous section. 



Fig. 5. — Field lines of the magnetic configuration associated with a 
relaxed equilibrium obtained from the potential configuration shown in 
Fig. 1 after twisting boundary motions have been applied up to / = t s and 
relaxation up to / = to', case t s = 400t a . Although the configuration is made 
of a continuum of field lines, only field lines belonging to three flux tubes 
are shown, as well as outer arcades, so that the effects of central strong shear 
and outer twist appear. [ See the electronic edition of the Journal for a color 
version of this figure.] 


6.1. Boundary Motions 

Our boundary motions are chosen to converge toward 
the inversion line {.v = 0} at the velocity 


= -B : (x,y,0.0)y = -£(x,y)y, (92) 

which practically vanishes, as required for consistency, on 
the boundary of the basis of our box. This velocity is associ- 
ated with functions / and g , which are both nonzero (they do 
satisfy d x f = —d y g\ but they can be computed only 
numerically. A cut of this flow along the axis {x = 0} (Fig. 
6) shows that it is continuous and has the typical property 
of the flows usually introduced in axisymmetric situations. 

The resulting evolution is still assumed to be ideal 
(r) = 0). As shown below, it is characterized by the existence 
of two qualitatively different phases: a quasi-static phase for 
to < t < t c and a disruption phase for t c < t. 

6.2. Existence of an Equilibrium Phase 

The simulations show the existence of a time t c > to such 
that the field evolves quasi-statically through a sequence of 
force-free equilibria for to < t < t c . Magnetic energy is not 
monotonically increasing from the beginning of this phase, 
as shown in Figure 3. For the reference case (corresponding 
to t s = 0) it is even monotonically decreasing. However, we 
get W(t c ) - W(to) > 0 for all the values of t s ± 0. It is inter- 
esting to note that the amount of magnetic energy injected 
in the configuration during this quiet phase is relatively 
large compared to the amount that has been stored during 
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Fig. 6. — (a) Vector plot and ( h ) cut and variations along the >’-axis (cut at {.v = 20}) of the boundary converging flow applied to the footpoints of the 
configuration for t > to. 


the global shearing-twisting process with / = 0. However, 
we still found W(t) < W„(t) all along. 

Magnetic shear (defined to be the angle between the inver- 
sion line and the transverse field) increases continuously, as 
seen in Figures 7 and 8. The effect of the converging boun- 
dary velocity field is to effectively make the field more 
aligned with the inversion line, but by keeping the topology 
to be that of a highly sheared arcade core. 

Magnetic helicity is almost constant during this phase 
[i.e., H(t) = //(/()) for to < t < t c ]. This result is fully consis- 
tent with equation (49), which applies here as equations 
(41), (44), and (45) are satisfied. The cases U 50 and U 400 
show dips in the magnetic energy evolution. We checked for 
U 50 that this dip (/ = 280 t a ) corresponds effectively to the 
existence of a neighboring equilibrium as shown in Figure 9, 
when the driver has been switched off with the linear ramp. 
This explains the still small rising of magnetic energy, which 
keeps a constant value, while kinetic energy keeps decreas- 
ing. However, we found that the reached equilibrium has 
still an arcade-like topology with no twisted flux rope. The 
case of U 400 , which will be discussed below, remains in 
equilibrium in the rising phase (for magnetic energy) that 
precedes the dip. 

6.3. Existence of a Disruption Phase 

As shown in Figure 3, there exists a t c at which the config- 
uration U 1 ' experiences a transition toward a dynamic evo- 
lution, with a nonnegligible amount of magnetic energy 
being dissipated [W(t) < W(t c )\. As shown in Figure 10, 
this transition is associated with a change in the topology of 
the magnetic field from an arcade-like form to a twisted 
ropelike form. This structure rises up, reconnecting strongly 
with the overlaying arcade lines, as well as below where new 
small arcades are observed to reform as shown in Figure 1 1 . 


Fig. 7.— Field lines showing the configuration obtained for t s = 200t a 
after a converging velocity field has been applied from to = 400r A to 
t = 450t a . [See the electronic edition of the Journal for a color version of this 
figure .} 
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Fig. 8.— Field lines showing the configuration obtained for t s = 200t a 
after a converging velocity field has been applied from / 0 = 400r A to 
t — 480t a (close to the nonequilibrium critical point). [See the electronic 
edition of the Journal for a color version of this figure.] 


Fig. 10. — Field lines showing the configuration obtained for t s — 200t a 
after a converging velocity field has been applied from /o - 400t a to 
/ = 498t a . The nonequilibrium critical point has been past, and magnetic 
reconnection leads to the formation of a three-part structure: twisted flux 
rope of about 27 t twist, small new arcades below, and overlaying arcade 
revealing a plasmoid-like magnetic structure. This structure is, however, far 
from equilibrium. [Scr the electronic edition of the Journal for a color version 
of this figure.] 



Fig. 9.— Evolution of the magnetic and kinetic energies (in units of W n \ 
see text) after the converging boundary velocity field has been switched off 
(with a linear ramp of / = 10r A ), in the dip (at t = 280t a ) of the energy 
curve of Fig. 3 for the case t s ~ 50t a . The existence of a neighboring equili- 
brium state is associated with the constant value reached by the magnetic 
energy while the total kinetic energy keeps decreasing toward zero. [See the 
electronic edition of the Journal for a color version of this figure.] 


A three-part magnetic structure is associated with this state: 
(1) a twisted flux rope running through (2) a global arcade 
and above (3) small loops. It is worth noticing that this 
reconnection going on (with transfer of twist) from the 
strongly sheared inner field lines toward the outer ones (less 
sheared) may even lead to a point where the center of the 
spots (corresponding to infinitely twisted lines) leads to field 
lines having a high twist. However, this state is not in equili- 
brium, as we obtain eventually v > p A in some region. 

As seen in Figure 12, the transverse magnetic field compo- 
nent exhibits a reversal in the sign of B v along the inversion 
line that is compatible with the appearance of the twisted 
flux rope, while the magnetic configuration remains arcade- 
like with high magnetic shear prior to the disruption. 

For all the cases considered in this study, we found that 
there is no neighboring equilibrium close to U l for some 
/ > t c . This point is checked by performing a relaxation pro- 
cedure during which the driving boundary velocity field is 
switched off. The configuration does not relax toward an 
equilibrium but still experiences a strong dissipation with 
acceleration of the plasma. Figure 1 3 shows the evolution of 
energies after the velocity has been switched off just after the 
critical point of disruption t s = 498ta. It is worth noticing 
that the configuration, although far from equilibrium, 
exhibits a twisted flux rope. 

The case t s = 400t a shows a first dip in the magnetic 
energy at / = 870ta. This dip does not correspond to the 
existence of a neighboring equilibrium, as shown in 
Figure 14, where the boundary driver has been switched off 
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Fig. 1 1 . — Field lines showing the configuration obtained for t s = 200t a 
after a converging velocity field has been applied from / 0 = 400t a to 
/ = 530t a . Reconnection goes on with the three-part structure conserved. 
One field line having a twist greater than lit between the arcade and the flux 
rope is shown. [See the electronic edition of the Journal for a color version of 
this figure.] 

with a linear ramp of 10r A . Magnetic reconnection associ- 
ated with a disruption occurs with a strong decrease in mag- 
netic energy, while kinetic energy rises strongly before 
slowly approaching a finite value. We have found that the 
magnetic topology during this disruption phase exhibits a 
twisted flux rope. It is worth noticing that unlike the other 
cases, the rate of decrease of magnetic energy becomes 
dominated by the injection rate associated with the converg- 
ing motions when they are present (as shown in Fig. 3), such 
that a second disruption phase occurs. However, during the 
second rising phase of magnetic energy, no equilibrium 
exists while the configuration shows a twisted flux rope far 
from equilibrium. 

It is worth noticing that, although for theoretical reasons 
we presented cases for which the velocity field is applied all 
along the evolution, this leads to unrealistic final configura- 
tions: at the end the two spots have been so much com- 
pressed against each other that they are almost reduced to 
two straight segments. This extreme limit may explain the 
small deviations of the magnetic helicity from keeping a 
constant value. It is interesting to remark that the existence 
of a critical time beyond which no relaxation toward a 
neighboring equilibrium exists implies that one does not 
need to approach this absolute limit. 

7. DISCUSSION 

In this paper we have addressed the problem of the evolu- 
tion of a three-dimensional magnetic field driven by con- 


verging motions imposed on the photosphere S, with the 
initial field Bo being a stable force-free field having a non- 
zero relative helicity. In a first step, we have considered a 
simple analytical model describing a situation of that type. 
In that model, the driving motions are of a quite peculiar 
type, as all the points of S suffer at any time / a homologous 
contraction toward the origin. The evolving field B(t) exhib- 
its, however, some interesting features: ( 1 ) B evolves quasi- 
statically through a sequence of force-free configurations 
that stay ideally stable and even unconditionally nonlinearly 
stable if Bo itself has this property, (2) the magnetic energy 
increases monotonically in time, and (3) the helicity keeps 
its initial value. This suggests that the system can suffer a 
disruption only if reconnection is allowed by a nonzero 
resistivity of the plasma. It is worth noticing that this simple 
analytical model may also be relevant to the divergence of a 
magnetic structure from the center of a supergranulation 
cell toward its borders, a characteristic picture of the quiet 
Sun. 

Next, we have reported the results of extensive numerical 
simulations. After having constructed a series of initial sta- 
ble force-free fields Bo = B(t 0 )y with |//(/ 0 ) | > 0 and still the 
initial B z on S, by deforming a given potential field in a two- 
step (twisting+relaxation) process, we have evolved these 
Bo by imposing on S motions converging toward their 
straight polarity inversion line. In all the cases we have con- 
sidered (either with small or large helicity), the evolution of 
the field exhibits two different phases. In the first one, the 
evolution is almost quasi-static. The magnetic topology 
remains arcade-like, with a shear along the inversion line 
increasing, magnetic energy being stored, and helicity keep- 
ing its initial value. The evolution in this phase is thus at 
least qualitatively quite similar to the one shown by the solu- 
tion of the analytical model introduced above. At some crit- 
ical stage, however, this quiet phase is stopped and the 
configuration experiences a transition to a second dynamic 
and strongly dissipative phase, during which reconnection 
leads to the formation of a twisted flux rope, however not in 
equilibrium. 

For defining the velocity field on the boundary (either 
during the twisting phase allowing to construct the initial 
field or during the converging phase), we have used a two- 
dimensional Helmholtz’s representation of B z v s (or equiva- 
lently of the tangential electric field) on 5, in which appear 
two scalar functions /and g. We have derived new formulae 
for the variation of the energy and the helicity in terms of 
these functions and checked in our calculations that the hel- 
icity injection rate was in good accord with that analytically 
derived. We have also checked that the magnetic energy 
was, at any time, smaller than the energy of the open field 
with the same flux distribution on the boundary. The con- 
jecture formulated by Aly ( 1 984) thus appears to be compat- 
ible with our calculations. 

Our three-dimensional results complete the two- 
dimensional ones reported earlier (Priest & Forbes 1990; 
Forbes & Priest 1995; Forbes 1991; Inhester et al. 1992). 
These authors have shown the existence of a catastrophic 
nonequilibrium transition, implying the ejection of a plas- 
moid, when the flux distribution evolves in some adequate 
way on the boundary. Their conclusions, however, were lim- 
ited by the existence of an unanchored flux rope and the 
unknown issue of nonequilibrium in three dimensions. This 
paper settles these important issues, showing that the system 
dynamically reconnects as soon as it loses equilibrium, 


No. 2, 2003 


CORONAL MASS EJECTION INITIATION. I. 


1085 






20.0 

col 


- 3 . 3029+00 — ^ - 5 . 3949+00 

(c) ( d ) 


Fig. 12— Evolution of the transverse magnetic field at {z = 0} for the case t s = 200t a at (a) t = t = 450t a , (/>)/ = / = 480t a , (c) t = / = 498t a , and (cl) 
t — 530 t a . The .v-axis O’-axis) runs along the direction named “ row ” (“ col ”). the electronic edition of the Journal for a color version of this figure. ] 


with no secondary intermediate nonequilibrium bifurcation 
being produced. 

Our results may be relevant to the problem of the initia- 
tion of CMEs, as they do show that a global disruption may 
occur in a magnetic structure with nonzero helicity contents 
when it is driven into an evolution by the converging 
motions that have been shown by some observations to be 
actually present on the photosphere. Let us comment briefly 
on two different important points, the origin of the helicity 
and the role of the flux rope, in these events. (1) As we have 
found that helicity keeps a constant value during the quasi- 
static phase of evolution, it needs to have been produced 
during a prior phase. In our simulations, helicity was 
obtained by twisting motions, but we do not claim that this 
process actually occurs in the corona; observational eviden- 
ces for adequate transverse photospheric velocities are still 
needed. It could as well be a result of processes taking place 
before the emergence of the structure as also estimated 
recently by Demoulin et al. (2002) and Nindos & Zhang 
(2002). We should also note that we cannot exclude that 


some amount of helicity be produced by the converging 
motions themselves, if they are less symmetric than the ones 
we have considered in our model. (2) Helical structures 
associated with prominences ejected as part of the CMEs 
are sometimes observed (but they have not been considered 
in this paper), and it is clear that twisted ropes are good can- 
didates for the support of cool material. It is still an open 
problem, however, whether a rope does exist prior to the 
disruption, thus possibly playing a role in its triggering. Pre- 
vious three-dimensional results had shown that both a 
sheared complex topology configuration of the multiarcade 
type (Antiochos et al. 1999) and a twisted flux rope (Amari 
et al. 2000) in a nonnecessarily bipolar configuration are 
candidates for the initiation of a CME. The results pre- 
sented in this paper complement these earlier ones, by show- 
ing an example of an evolving bipolar configuration 
suffering a major disruption, but without the presence of a 
twisted flux rope in equilibrium. A rope is created, but only 
as a result of reconnection during the global disruption, and 
it is then part of a highly nonequilibrium process. The 
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Fig. 13. — Evolution of the magnetic and kinetic energies (in units of W ff ; 
see text) after the converging boundary velocity field has been switched off, 
right after the nonequilibrium point, at t = 500t a seen on the energy curve 
of Fig. 3 for the case t x = 200rx. There is no relaxation toward a neighbor- 
ing equilibrium even at this stage. [See the electronic edition of the Journal 
for a color version of this figure.] 
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Fig. 14. — Evolution of the magnetic and kinetic energies (in units of W-\ 
see text) after the converging boundary velocity field has been switched off 
(with a linear ramp of t = 10r A ), in the dip (at t = 870t a ) of the energy 
curve of Fig. 3 for the case t s = 400t a . No neighboring equilibrium state 
exists, and a disruption occurs. [See the electronic edition of the Journal for a 
color version of this figure .] 


results presented here are actually comparable to those 
obtained for axisymmetric sheared arcades, in which there 
is a high shear all along the inversion line. Such a shear on a 
large coherent length scale, impossible to reach in three 
dimensions by pure vortex motions, has been obtained here 
thanks to converging motions. 


We acknowledge support from NASA's Sun-Earth 
Connection Theory Program, NASA's STEREO/ 
SECCHI Consortium, and Centre National d'Etudes 
Spatiales, which also supported Dr. Amari's visits to 
SAIC in San Diego. The numerical simulations per- 
formed in this paper have been done on the NEC SX5 
supercomputer of the IDRIS of the Centre National de 
la Recherche Scientifique. 
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Introduction 

Simulating CMEs on active-region length scales is challenging 

We have been developing a 3D MHD model for this task 

• Non-uniform meshes 

• 3D finite differences in spherical (r,0,0) coordinates 

• Implicit and semi-implicit time differencing 

• Comprehensive physics model including the solar wind and 
energy transport (radiation, parallel thermal conduction, heating, 
and Alfven waves) 



Evolution of an Active Region in 1996 

Kitt Peak Magnetogram Big Bear H-a Yohkoh SXT 

Aug 29. 1996 


Sep 26, 1996 Sep 25, 1996 Sep 25, 1996 



Figure 4. Magnctograms, l ull-disk H -a images, and Yohkoh soft X-ray images of an active region in laic August 
and late September 1996. The two dales arc separated by approximately one solar rotation. Note that the H-« 
images show the presence of a filament. Note also that the magnetic field in the active region has diffused and 
dispersed from August 29 until September 26. The SXT image shows an S-shapc in the active region. 
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Figure I. Modeling the magnetic and thermal structure of an acli\c region on August 29, 1996. A Kitt Peak 
magnetogram is used to specify the normal component of the magnetic field. A tw ist is applied to the field, and the 
steady state is calculated for a given coronal heating distribution. The temperate and dcnsit\ structure shows that the 
transition region height varies in different parts of the active region. 
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Example 1: 3D Streamer 

• A typical coarse run takes a few hours on 8-64 processors 

• 71x51x51 mesh on 32 processors on the IBM/SP3 took 2 hours 

• Nonuniform meshes with: 

Armax/Armin ~ 70, A0max/A0rnin ~ 5, A0max/A0min ~ 10 

• Dipole magnetic field aligned with the solar rotation axis, together 
with an added sub-surface dipole to model an active region 

• Solar wind with a polytropic energy equation with y = 1 .05 and a 
temperature of 1 .8 x 10 6 K and a density of 10 8 cm' 3 in the lower 
corona 






Example 2: Shearing a Model Active Region 

• Dipole magnetic field aligned with the solar rotation axis, together 
with an added sub-surface dipole to model an active region 

• Start with a potential field 

• Apply a shearing flow that preserves B r : 

Vshear(0,0) = Vj_ X IpY 

where ip(d,(p) = B, 2 (r = Ro,d,(f)) 

• This flow energizes the magnetic field and builds free magnetic 
energy 

• This is in preparation for CME initiation (by continued shear or by 
flux cancellation) 


Localized Shearing of an Active Region 


Radial Magnetic Field Shear Flow 




Shearing an Active Region in a Localized Manner 



Potential Field 



Sheared Field 


May 12, 1997 Event 



May 1 1 , 04:05UT 


Figure 5. A sequence of MD1 magnetograms showing convergence of opposite polarity magnetic near the neutral 
line flux prior to the launch of the CME. 


September 12, 2000 Event 


SOHO/MDI Magnetograms 



Figure 7. A sequence of MDI magnetograms for the September 12, 2000 CME. The complex nature of the 
magnetic field is apparent, as the Sun is approaching solar maximum at this time. The magnetic evolution shows 
some c\ idencc of flux cancellation and dispersal of the magnetic Held in the active region prior to eruption. There is 
also significant surrounding flux which could lead to rearrangement of the outer field lines. Either the flux 
cancellation model or the breakout model may be consistent with this event. 



A Model Active Region (7986 of August, 1996) 
(interpolated and smoothened) 

1. A force-free field, generated by shearing the base with a flux- 
preserving flow, is shown below (top view). 

2. The torsion parameter a of the structure is approximately the 
same as the one obtained by Madrini et al. (Geofisica 
International 2000, 39 ,73). 



Important Assumptions and Consequences 


1. Magnetic field is sufficiently strong so that the thermal 
conductivity is finite only along the field. 

2. Plasma (3 is sufficiently small so that the field is relatively rigid 
and the plasma dynamics is essentially one dimensional. The 
magnetic structure behaves like an ensemble of static 
individual field lines. 

3. Three models for the plasma heating rate have been 
investigated: 

(A) heating - B 1 2 3 4 , stochastic energy buildup, current layer 

dissipation 

(B) heating - Vaifven, magnetic reconnection 

(C) heating - density 2 , a phenomenological model using EIT 

data fit of isothermal loops (Neupert et 
al. 1998, SP. 183, 305). 

4. For each model, we solve the system of dynamic-energy 
equations to reach a steady state in order to obtain the 3D 
profiles of temperature and density. A deep chromosphere in 
the computation serves as a mass reservoir. 
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SIMULATION OF CMES ORIGINATING IN 
ACTIVE REGIONS 

Z. Mikic. J. A. Linker, P. Riley, and R. Lioncllo 

Science Applications International Corporation, San Diego, California, USA 

Previously we have addressed Ihc initiation of CMEs by large-scale changes in the 
solar phoiosphcric magnetic flux. These CMEs had a global nature, since their size 
was comparable to the solar radius. These early investigations were primarily focused 
on the theoretical aspects of CME initiation, and the models were thus rather idealized. 
We will describe recent advances in our computational models that have extended our 
capability to study localized CMEs, such as those that might originate from an active 
region. This capability will allow us to analyze actual CME events, and to compare 
our results in detail with CME observations. 


Research supported by NASA and NSF. 



MHD Equations 

(Improved Energy Equation Model) 
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= J x B - V/7 - Wp w + pg + V-(vpVv) 
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X + V-(pv) = (y - l)(-pV-v - V-q - n e n p Q(T) + H) 


y =5/3 

A /V 

q = -^,bb-Vr (Close to the Sun, r < 10 R s ) 

A /V 

q = 2an e T bb*v/(y - 1) (Far from the Sun, r > 10 R s ) 
+ WKB equations for Alfven wave pressure p w evolution 


Code Properties 

Written in FORTRAN 90 

Designed to run on massively parallel computers using MPI 

• Linux & Beowulf (lf95, pgf90, Intel Fortran) 

• Mac (Absoft) 

• IBM/SP3 (xlf) 

Mesh decomposition among processors in 3D 

We do not use FFTs in 0, so that the mesh can be nonuniform 

Dynamic allocation allows mesh size and number of processors to 
be selected at run time 

Restart capability using HDF files (for long runs) 

Essentially a complete re-write of our spherical 3D MHD Cray 
code MAS 



Conclusion 

The new 3D MHD code runs on parallel computers and is expected 
to scale well 

The code can be used to perform high-resolution runs (300 3 ) on 
massively parallel computers 

We have developed the techniques to localize mesh resolution in 
particular regions (albeit with structured meshes) 

This capability will allow us to simulate a compact CME 
originating in an active region 

With this capability we can study specific events, allowing detailed 
comparison with CME observations 

Three possible events that may be studied in the future are a solar 
minimum event (May 12, 1997), an event during the rising phase of 
the solar cycle (May 1998), and a solar maximum event 
(September 12, 2000) 
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Abstract 

The thermal structure of an active region depends on the mechanism that heats 
the coronal plasma. A number of coronal heating mechanisms have been 
proposed over the years. They have different parametric dependence on the 
magnetic field, plasma density, and possibly other variables. Different 
mechanisms result in different thermal structures, and therefore, different EUV 
and soft X-ray emissions from an active region. Hence, the comparison between 
the computed emissions based on these models and the observed emissions will 
help to discover the parametric dependence of the actual heating mechanism and 
put some restrictions on the theoretical models. We have developed a 3D thermo- 
magnetohydrodynamic code to compute the thermal structure of an active region. 
The emissions resulted from various heating models will be compared with the 
images obtained from SOHO and Yohkoh. 


Objectives of Study 

1. To construct the 3-D thermal structure of an active region 

2. To understand how plasma heating affects the thermal structure 

3. To compute the EUV and soft X-ray emissions from the structure 

4. To demonstrate that the technique we have developed can be used as a 
possible diagnostics for plasma heating in the corona 



Preliminary Results + Discussion 


1. Different heating models lead to substantially different thermal 
structures. In model A, the hot region is confined closely to the 
neighborhood of the magnetic poles, as the heating rate (-B 2 ) decreases 
rapidly with distance. It reaches only to a moderate altitude ( - 0.06 
solar radius). In model B, the hot region reaches to a higher altitude 
and spans a larger volume because the heating rate (- VAifven - B p 1/2 ) 
decreases more slowly compared to model A. In model C, the hot region 
is confined to a tiny region above the neutral line because of the rapidly 
decreasing density above the transition. 

2. The density profiles are generally consistent with the temperature 
profiles at equilibrium according to the gravitational scale height. 

While the temperature profiles appear to be smooth, the density 
profiles above the magnetic poles show fine structures. The fractional 
variation in temperature is smaller than the fractional variation in 
density. The fine structures are plasma loops whose field lines are 
longer than their immediate neighbors, a probable condition in a non- 
potential field. The longer field line reduces the effectiveness of thermal 
conduction, resulting in a slightly higher temperature and a longer 
gravitational scale height at equilibrium. 

3. Thin, EUV emitting loops appear naturally due to the elevated density 
and temperature of these plasma loops. The location and size of these 


thin loops are correlated with the high-density region in each model. In 
model A, the apex is at -0.06 solar radius. It is -0.09 in model B, and 
ranges from 0.01 to 0.03 in model C. In addition to the thin loops, there 
are some diffuse EUV loops that reach a higher altitude. 

4. The computed soft X-ray images show somewhat diffused features as 
expected. Some loop-like structures are visible. 

5. Limitations of this study — there are many - in this first attempt to 
obtain a 3D thermal structure of an active region: 

A. Vector magnetograms are not available for this case. Therefore, 
there are uncertainties about the geometry of the field lines, 
which critically determine the direction of heat flow. Future 
technology should provide more reliable information about the 
field so that the computed thermal structure is more realistic. 

B. The magnetogram has been highly smoothened and interpolated 
in the simulation so that it can be handled in a 127x89x127 mesh. 
The fine structures in the field are lost. Therefore, it is not 
possible to reproduce the multiple fine loops as seen in the actual 
SOHO/EIT and Yohkoh/SXT images. The original, unprocessed 
magnetogram would have required more grid points than that 
can be handled practically. 

C. The sharp density-temperature gradients in the transition region 
require very fine vertical mesh to resolve, making the 
computation very time consuming. As a result, approximations 
must be made to make the computation practical and within our 


computing resources. We have attempted to minimize the 
numerical errors, but it can still be improved. The code is 
being ported to massively parallel environment, and faster 
hardwares will greatly improve the accuracy of the results. 

The computed images from the three models cannot consistently 
match the actual ones across all wavelengths. 

This is not completely surprising because of the uncertainties in 
the input and the approximations we have to make in this first 
attempt. In addition, other plausible models have not been 
considered. Nevertheless, we have developed the framework and 
the tools needed for future works that will address these issues. 
Stay tuned. 


I 


Emissions from the Computed 
Thermal Structures 

1. EUV emissions at 171, 195 and 284A are computed by using the 
SolarSoft routines and the CHIANTI Atomic Database. The 
instrument response of EIT has also been taken into account. In 
each case, the emission is integrated along the line of sight. For 
each heating model and each line, we have computed the 
expected emission from two different perspectives, a top view as 
if the structure is at the center of the disk and a side view as if it 
is at the limb. 

2. Soft X-ray emission is also computed from the same perspectives 
and is integrated along the line of sight. The Yohkoh/SXT 
instrument response has also been taken into account. 

3. In the panels below, we plot the temperature and density 
profiles at cross section A for each of the heating model to show 
the details of their vertical structure above the active region. To 
their right are the EUV and soft X-ray emissions from two 
perspectives. The familiar SOHO/EIT color scheme is used. 
(Blue=171A, Green=195A, Yellow=284A). 
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